Velocity distributions in dilute granular systems 
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Motivated by recent experiments reporting non-Gaussian velocity distributions in driven dilute 
granular materials, we study by numerical simulation the properties of inelastic gases as functions of 
the coefficient of restitution r\ and concentration cj> with various heating mechanisms. We show that 
there are marked, qualitative differences in the behavior for uniform heating (as is frequently assumed 
theoretically) and for particle systems driven at the boundaries of the container (as is frequently 
done in experiments). In general, we find Gaussian velocity distributions for uniform heating and 
non-Gaussian velocity distributions for boundary heating. Furthermore, we demonstrate that the 
form of the observed velocity distribution is governed primarily by the coefficient of restitution n 
and q — Nh/Nc, the ratio between the average number of heatings and the average number of 
collisions in the gas. The differences in distributions we find between uniform and boundary heating 
can then be understood as different limits of q, for g > 1 and q < 1 respectively. Moreover, we 
demonstrate that very similar behavior is found for a simple model of a gas of inelastic particles 
with no spatial degrees of freedom. 

PACS numbers: 81.05.Rm, 05.20.Dd, 83.10.Pp 



I. INTRODUCTION 



Granular materials consisting of macroscopic particles 
or grains can exhibit behavior reminiscent of conventional 
phases of matter. Sand, for instance, can flow like a liquid 
under some conditions. Dilute granular systems, or gases, 
have been extensively studied both experimentally and 
theoretically, in large part as simple model systems ex- 
hibiting noncquilibrium and dissipative behavior. These 
systems are intrinsically dissipative and out of equilib- 
rium, even though it is tempting to apply such equi- 
librium notions as temperature. Since the collisions in 
such a gas are inelastic, it is necessary to supply energy, 
or to drive them in order to maintain a gas-like steady 
state. Otherwise, inelastic collapse can occur, in which 
all motion ceases after only a finite time 0,0- In princi- 
ple, it is possible to heat or drive the system uniformly 
throughout the container (uniform heating), as has been 
done in simulations 0, Q, and as is assumed in analytic 
theories 5). In experiments, however, one usually drives a 
granular gas by shaking or vibrating the walls of the con- 
tainer. Such boundary heating means that the ener gy i s 
inserted in a spatially inhomogeneous wav|rj|^.l8ll^llC|. 
As a consequence, the gas will develop a gradient in den- 
sity and mean kinetic energy. 

Even for uniform heating, however, significant devia- 
tions from equilibrium gases, e.g., in density correlations, 
are observed 3j. It is often assumed that in the bulk, 
where the boundaries are far away and spatial gradients 
are small on the scale of the mean free path of the par- 
ticles, both heating methods should give the same be- 
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havior. Here, we show that there are striking differences 
in the behavior of granular gases heated uniformly or 
through the boundary. 

There has been much interest recently in another as- 
pect of these non-equilibrium gases: namely, the velocity 
distributions. This is in part because they deviate from 
the Gaussian distributions that one would expect if the 
collisions were elastic. It has been suggested that the ve- 
locity distributions can be described by a sort-of stretched 
Gaussian, of the form P(v) = Cexp[— (3(v/a) a ), where 
a = (i> 2 )2 is sometimes called the granular temperature. 
Rouyer and Menon have reported such a distribution 
with exponent a = 1.5 over the whole observed range 
of velocities, which was unaffected by changes in am- 
plitude and frequency of driving. This observation was 
particularly intriguing, given the theoretical predictions 
obtained before by Van Noije and Ernst They devel- 
oped a kinetic theory that predicts a high-velocity tail 
described by a distribution with an exponent a = |. 

Unfortunately, results for velocity distributions are 
never unambiguous. Both in simulation and experi- 
ment, different setups and driving mechanisms usually 
give different behavior of the velocity distribution. For a 
setup where particles on a horizontal plate were driven 
in the vertical direction, Olafsen and Urbach [llj found 
a crossover from exponential to Gaussian distributions 
as the amplitude of the driving was increased. The re- 
sult of Rouyer and Menon pfl was obtained for a different 
configuration where particles were confined between two 
vertical plates and driven in the vertical direction. Al- 
though the exponent a = 1.5 they find is reminiscent 
of theoretical results obtained by Van Noije and Ernst 
both results are actually inconsistent in that the ex- 
ponent a — | can only describe the high-velocity tail 
of the distribution. Such a high-velocity tail was actu- 
ally obtained by Moon et al. 0. Using a simulation 
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with a variation on uniform heating they find a velocity 
distribution that has an exponent of a = 2.0 for low ve- 
locities, but crosses over to an exponent of a = 1.5 for 
the high-velocity tail. The results of Rouyer and Menon 
are also limited in the sense that they were obtained for 
almost elastic particles with rj = 0.93. Blair and Ku- 
drolli 8] use a different setup where particles move along 
an inclined plane. Friction with the plane during col- 
lisions reduces the coefficient of restitution to 77 w 0.5. 
They find the distribution with exponent a = 1.5 only 
in the very dilute case. Otherwise, the distributions de- 
viate strongly from both Gaussian and the distribution 
obtained by Rouyer and Menon. It is interesting to notice 
that for the denser case Blair and Kudrolli find a distri- 
bution with a crossover much like the type observed by 
Moon et al. 

At present there is no agreement on what the veloc- 
ity distribution of a granular gas looks like exactly nor is 
it clear what causes the velocity distributions to deviate 
from a Gaussian. Puglisi et al. have suggested that 
the deviations are caused by the spatial correlations in 
the gas. They propose that for a region of uniform den- 
sity the velocity distribution of the gas actually is Gaus- 
sian, with a density dependent width. The spatial cor- 
relations cause density fluctuations and they claim that 
the non-Gaussian distributions arise as an average over 
the velocity distributions over these regions of different 
density. It is true that an average over Gaussian distribu- 
tions with different widths in general yields a stretched 
Gaussian, but experiments performed by Olafsen and Ur- 
bach ^l| showed that this does not happen in granu- 
lar gases. They find that even within small windows of 
uniform density the velocity distributions remained non- 
Gaussian. Here, we show that non-Gaussian behavior 
arises even in a simple model with no spatial degrees of 
freedom. 

We study behavior of the velocity distributions of 
the granular g function of </>, the area fraction, 

and 77, the coefficient of restitution. Specifically, we 
consider the effect on the velocity distribution of driving 
the gas by heating uniformly, as is assumed in theory 
and many prior simulations, and by heating through a 
boundary, as is done in most experiments. We will show 
that there exists clear qualitative difference between the 
velocity distributions for uniform and boundary heating. 
Furthermore, we will show that there is no evidence for a 
universal velocity distribution with a constant exponent 
a = 1.5. We demonstrate instead that a family of dis- 
tributions showing apparent exponents covering a wide 
range of values a < 2 is expected, depending on both 
material and experimental conditions. Furthermore, we 
show that the velocity distribution is governed primarily 
by the relative importance of collisions to heating, i.e., 
the way in which energy flows through the system of 
particles. Specifically, we introduce a new parameter 
q = Nh/Nc, which measures the ratio between mean 
numbers of heating events Nh and mean number of 
collisions Nc experienced by a typical particle. These 



theoretical observations can explain both the observed 
non-Gaussian behavior, as well as the ambiguities in the 
experimental and theoretical literature on dissipative 
gases to date. We also demonstrate that the behavior of 
the velocity distributions can be captured quantitatively 
by a simple model that takes only 77 and q into account, 
with no spatial degrees of freedom. 



II. NUMERICAL SIMULATION 

We use an event-driven algorithm to simulate N par- 
ticles of radius r moving in a two-dimensional box. Par- 
ticles gain energy by heating and lose energy through 
inelastic collisions. When two particles i and j collide 
their final velocities depend on their initial velocities in 
the following way: 

v i = v i — ( v * ' r ij ~ Vj • ry)ry , (1) 

where < 77 < 1 is the coefficient of restitution and f ^ is 
the unit vector connecting the centers of particles i and 
j- 

For uniform heating we adapted an one-dimensional 
algorithm described in [j| . When heating uniformly, each 
individual particle is heated by adding a random amount 
to the velocity of each particle during a time step At: 

Vi (t + At) = Vi + VhAtf(t) , (2) 

where f(t) is a vector whose components are uniformly 
distributed between — ^ and ^ and h is proportional to 
the heating rate. After heating the system is transferred 
to the center-of-mass frame. Particles move in a box with 
sides L = 2. We use periodic boundary conditions to sim- 
ulate bulk behavior. The time step At is chosen in such 
a way that on average the number of collisions per time 
step is less than one. It should be noted that this heat- 
ing mechanism is significantly different from the s pat ially 
homogeneous heating used in some experiments |ll|. In 
the experiment all particles feel the same forcing, so the 
motion of the neighboring particles is strongly correlated 
in space and time. In uniform heating however, all parti- 
cles are independently driven by a stochastic source and 
as a consequence correlations are very weak. 

When heating through the boundaries, particles gain 
velocity upon collision with the boundary. For simplicity, 
we assume that the collision with the boundary is elas- 
tic. In that collision occurs by reflecting vj_, the 
component of the velocity perpendicular to the boundary. 
Heating occurs by adding a random amount of velocity 
to v_i_. Then after collision with the boundary one has: 

= v - 2v_l + Vhf (i). (3) 

Particles move in a circular box of radius R = 1. A 
symmetrical container has the advantage that it allows us 
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to examine density and granular temperature gradients 
along a single coordinate r, the distance from the center 
of the box, as in the one-dimensional case Q . 

We start the simulation by distributing the particles 
uniformly over the box. When using boundary heating, 
we give each particle a small, uniformly distributed ve- 
locity to enable particles to reach the boundary. Then 
particles are heated and we allow the system to reach a 
steady state before taking data. For both uniform heat- 
ing and boundary heating, data is taken periodically ev- 
ery time step At. For uniform heating, data is taken 
when the particles are heated, so At equals the time be- 
tween heating events. 



III. SIMULATION RESULTS: CLUSTERING 

One of the first striking differences between uniform 
heating and boundary heating is clustering. When heat- 
ing through the boundary, a stable liquid-like cluster sur- 
rounded by a hot gaseous state will form for low coeffi- 
cients of restitution 77 or high area fraction (j). This occurs 
as particles are compressed in the center of the box by 
the pressure of particles moving in from the boundary. As 
the cluster grows in size, it can no longer be destroyed 
by the impact of high velocity particles and the cluster 
remains stable. A typical example is shown in Fig.^ We 
do not find these clusters when heating uniformly. This 
is because all particles are heated all the time, which 
prevents the collapse to a cluster. In experiments that 
find clusters in homogeneously driven systems 11] this is 
different, because here the driven motion of neighboring 
particles is highly correlated. The formation of a cluster 
spoils measurement of the velocity distribution because 
the system is no longer in a pure gas-like state. To avoid 
values of cj> and r\ corresponding to the formation of clus- 
ters in our simulation, we constructed a phase diagram. 
We did this by counting for every particle the number 
N$ r of neighbors with their center within a distance 6r 
from that particle. When the gas is in a hexagonal close 
packed state N$ r = 32. We obtained the distribution 
P(Ne r ) for different values of <f> and rj. An example for 
N = 350 and cf> = 0.1 is shown in figure |2 

For r] = 0.9 the distribution corresponds to a state with 
the particles uniformly distributed over the box and the 
peak of the distribution at the mean value N& r = 3.6. 
For r\ — 0.7 the distribution becomes bimodal, with a 
broad peak at high N 6r corresponding to the densely- 
packed cluster and a peak at Ne r = 1 corresponding to 
the surrounding dilute gas. The distribution shows a 
continuous transition for rj in between, which makes it 
hard to pinpoint an exact value of rj for which the gas 
enters the clustered state. Still, by looking at the shape 
of the distributions, it can be argued that the transition 
occurs somewhere between 77 — 0.75 and 77 — 0.85. This 
was repeated for different values of 4>, which allowed us 
to determine a sort of phase diagram. Specifically, we 
determined the limit of a pure gas-like phase, and all 



results presented below were obtained in this state. 



IV. SIMULATION RESULTS: VELOCITY 
DISTRIBUTIONS 

The velocity distributions P(v x ) for uniform heating 
are shown in Figs. El and ^| The velocity component 
v x is scaled by a x — (w^) 5 and the maximum of the 
distribution P{v x /a x ) is scaled to be unity. For a broad 
range of the parameters <fi and 77 the velocity distributions 
are very close to Gaussian. For 77 = 0.8 the velocity 
distributions can be fitted by a distribution with a = 
2.0. This decreases only slightly for 77 = 0.1, which can 
be fitted by a distribution with a — 1.9. Values of a 
are found to be independent of <fi. These exponents are 
constant over the entire observed range of velocities and 
we see no high- velocity tails with a — 1.5 for the range 
of <fi and 77 we examined. 

For boundary heating the gas develops a gradient in 
both density and mean kinetic energy as shown in figures 
and El Ideally, we want to measure velocity distribu- 
tions in a region where the gradient is small. To this end 
we divided the box in five rings of width 0.2. These rings 
are indicated in figures [5] and Only for values of </> and 
77 close to the clustering state, does the density within a 
ring vary by more than 10%. The velocity distributions 
P(v x ) for particles within the different rings are shown 
in figure [3 Figure Ha) shows P(v x ) with the velocity 
component v x scaled by a x = (w^) 5 and the maximum 
of the distribution scaled to be unity. All distributions 
seem to have the same shape. This is remarkable, since 




FIG. 1: Snapshot of a clustered state for N = 350, cj> = 0.05 
and 77 = 0.6. The circles indicate the current positions of the 
particles, while the lines indicate the direction and magnitude 
of their velocity. 
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FIG. 2: Number of neighbors within a distance 6r of a given 
particle for N = 350, <f> = 0.1 and 77 = 0.7(o), 0.75(D), 0.8(o), 
0.85(A) and 0.9(<). On average N 6r = 3.6 for <j> = 0.1. 
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FIG. 3: (a) P(v x /a x ). (b) -bn{-]n[P(v x /a x )]} versus 
\n(v x /a x ). Data for both figures is taken for TV = 350 and 
for = 0.02 and ij = 0.8(o), 0.6(D), 0.4(o), 0.2(A), 0.1(<). 
A Gaussian is shown as a solid line with slope —2 and the 
distribution obtained by Rouyer and Menon is shown as the 
dashed line with slope —1.5. 



naively one would expect the dynamics of the particles 
to be different for the different rings. In the center the 
density is high and particles have many collisions, causing 
strong spatial correlations. Along the boundary the gas is 
dilute. Here, many particles have experienced recent col- 
lisions with the boundary, and spatial correlations should 
be much weaker. It is has been suggested that these cor- 
relations cause the velocity distributions to deviate from 
Gaussian. But contrary to the expectations, within the 
container the shape of the distributions seems to be de- 
termined by the granular temperature a x only. This is 
a feature that is observed for all values of 4> and r), even 



close to the cluster state. 

Figure[7Jb) shows the behavior of the exponent a. This 
behavior is very different from the case of uniform heat- 
ing. For uniform heating a has the same value over the 
entire observed range of velocities. For boundary heat- 
ing, on the other hand, a has a constant value a± over 
the low-velocity range but crosses over to different value 
a2 when above a critical velocity v c . For the two inner- 
most rings a\ = 1.7 and «2 = 1-0. For the third ring 
ot\ = 1.8 and the outer ring has a = 1.5. For the outer 
two rings the velocity distributions actually undergo two 
crossovers. First it crosses over to a 2 ~ 1.0 and then 
the exponent increases again to ~ 1.5. In figure we 
show the effect of a change in cf> and 77 on the shape of 
the velocity distributions. Here we focus on the velocity 
distribution as measured in the ring with 0.4 < r < 0.6. 
This has the advantage of good statistics, but for values 
of 4> and 77 close to a cluster, we might see effects due to 
the density gradient in the gas. As shown in figure Ufa) 
the exponent ot\ = 1.8 except for 77 = 0.4, where ot\ = 1.6. 
For 77 = 0.9 there is no crossover in the observed range 
of velocities. In the other distributions one does observe 
a crossover and the point where it occurs shifts down to 
lower velocities as 77 is decreased. It is clear that the 
distribution for velocities above the crossover cannot be 
described by a single exponent. For low enough 77, the 
distribution seems to approach a constant exponent for 
high velocities. This exponent decreases from a 2 = 1.3 
for 77 = 0.7 to a 2 = 1.0 for 77 = 0.4. In figure G^b) 
a-i = 1.8 for (j) = 0.01 and <j> = 0.02, a x = 1.9 for = 0.03 
and cti = 1.7 for <f> = 0.05. For every <fi does one observe a 
crossover and the velocity at which the crossover occurs 
shifts only a bit as <f> is varied. The distributions ap- 
proach a constant exponent for high velocities. This ex- 
ponent goes down from a.% = 1.5 for (j> = 0.01 to «2 = 1-0 
for cf> — 0.05. In general, the deviations from Gaussian 
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FIG. 4: (a) P(v x /a x ). (b) -bn{-]n[P(v x /a x )]} versus 
\n(v x j <j x ) . The dashed lines have slope —2 and —1.5. Data 
for both figures is taken for iV = 350 and for 77 = 0.2 and 
<t> = 0.1(o), 0.05(D), 0.02(o), 0.01(A). 
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FIG. 5: (a) The average number density p as a function of 
distance r to the center of the box. Data taken for N = 350, 
cj> = 0.02 and 77 = 0.9(o), 0.8(D), 0.7(o), 0.6(A), 0.5(«) and 
0.4(t>) . (b) The mean kinetic energy (v 2 ) per particle as a 
function of r, for the same values of cj> an d V- The dashed 
lines indicate the concentric rings, within which the velocity 
distributions were separately calculated. 



become more pronounced as dissipation increases, i.e. as 
eft increases or as rj decreases. 

To test whether the velocity distributions we find here 
are only observed for this specific driving mechanism of 
heating through a circular boundary, we constructed dif- 
ferent systems that drive through boundaries in a differ- 
ent way. For instance, we constructed a box with periodic 
boundary conditions that includes a small circular region 
around the center. Within this circular region particles 
are uniformly driven but outside of the region they are 
not heated at all. For particles within the circular region 
we observe velocity distributions that are Gaussian. On 
the other hand, for particles outside of the circular region 
we observe the same non-Gaussian velocity distributions 
as seen in the case of a circular boundary. 

For uniform heating the distribution of velocities of 
particles that are heated is Gaussian, because the veloc- 
ities undergo a random walk. This is not the case for 
the particles heated at the boundary. To see if velocity 
distributions are sensitive to the precise way of heating 
at the boundary, we changed our heating algorithm so 
that, when a particle hits a boundary, it's new velocity is 
drawn from a Gaussian distribution. This has a modest 
effect on the far end of the high- velocity tail, but leaves all 
major differences between uniform and boundary heating 
intact. 

Finally, we studied the behavior of the velocity dis- 
tribution for different particle numbers N. Figure [Sfa) 
shows that, as N increases, the velocity distributions fall 
off more rapidly, but the high-velocity tails grow longer. 
This is shown more clearly in figure |5Jb). For the high- 



velocity tail we find a = 1.7 for N = 100 decreasing to 
a = 0.7 for N — 1000. The crossover shifts to lower ve- 
locity and becomes sharper as N is increased. As we will 
demonstrate in more detail below, for boundary heating 
there is no thermodynamic limit. Instead of approaching 
a limiting velocity distribution as N is increased, we find 
that the shape of the velocity distribution depends not 
only on 77 and 0, but also on N. 

For their experiments Rouyer and Menon used N par- 
ticles with 77 w 0.9, where 100 < N < 500 and 0.05 < 
eft < 0.25 0- In figure ITUI we plotted the velocity distri- 
bution for T) — 0.9, eft = 0.05 and several values of N. The 
solid black lines indicated the fit with a = 1.52 as made 
in |2J. This line clearly coincides with the high- velocity 
tail of the velocity distribution found by simulation. This 
suggests that it is possible that instead of a universal dis- 
tribution with a = 1.5, they only observed a part of a 
more complex velocity distribution, with two exponents 
and a crossover. For higher eft the high- velocity tails have 
exponents that are smaller than the a = 1.5 observed for 
eft = 0.05. 

The main difference between uniform and boundary 
heating is that in the first case heating takes place homo- 
geneously throughout the box, whereas in the latter case 
energy is injected inhomogeneously at the boundaries. 
This is not the direct cause for the difference in velocity 
distributions. This is demonstrated by the occurrence of 
non-Gaussian velocity distributions in a system studied 
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FIG. 6: (a) The average number density p as a function of 
distance r to the center of the box. Data taken for N = 
350, T) = 0.9 and cj> = 0.1(o), 0.05(D), 0.02(o) and 0.01(A). 
(b) The mean kinetic energy (v 2 ) per particle as a function 
of r, for the same values of eft and rj. Note that even for 
the dilute case eft = 0.01 the mean kinetic energy profile is 
not constant, but drops at the boundary of the box. The 
profile only becomes constant after a certain distance into the 
container, that corresponds to the mean free path of particles 
leaving the boundary. 
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FIG. 7: Velocity distributions calculated separately within 
the concentric rings shown in Figs. and |S| i.e., for < r < 
0.2(o), 0.2 < r < 0.4(D), 0.4 < r < 0.6(o), 0.6 < r < 0.8(A) 
and 0.8 < r < 1(<), where r is distance to the center. Data 
were taken for N = 350, cj> = 0.05 and rj = 0.8. (a) P(v/a x ). 
(b) — ln{— hx[P(v x /a x ))} versus ln(v x /<r x ). 
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FIG. 8: (a) -ln{-ln[P(v x /cr x )]} versus ln(v x /a x ) for N = 
350, <t> = 0.02 and 77 = 0.9(o), 0.8(D), 0.7(o), 0.6(A), 0.5(<) 
and 0.4(v). (b) — ln{— la[P(v x /a x )]} versus ln(v x /a x ) for 
N = 350, T) = 0.7 and <j> = 0.01(o), 0.02(D), 0.03(o) and 
0.05(A). 

by Moon et. al. 4J. They study a uniformly driven gas, 
but nevertheless find velocity distributions with a clear 
crossover. The main difference is the way they implement 
uniform heating. Every time At, they select at random 
two particles and add to these particles a random but 
opposite velocity to conserve the total momentum. On 
average heating is spatially homogeneous and in the limit 
of very small At the behavior of both their and our uni- 
form heating algorithm is the same. For small enough At 
the number of collisions between heating events is smaller 
than one and the heating dominates over the dissipative 
behavior of the gas. When At becomes larger, many col- 
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FIG. 9: (a) P{v x /a x ) for different values of N. (b) 
— ln{—\n[P(v x /a x )]} versus ln(v x /a x ). Data is taken for 
cj> = 0.05, v = 0.8 and TV = 100(o), 200(D), 500(o), 700(A) 
and 1000(<). 



lisions might happen between two heating events. Apper- 
ently this changes the velocity distributions. In figure HT1 
we show the velocity distribution for N = 350, <fi = 0.02 
and r\ = 0.4. The gas is heated using the two-point heat- 
ing algorithm described above, varying the time between 
heatings, At. For At = 0.01 the distribution has a ex- 
ponent a = 1.7 that is constant approximately over the 
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FIG. 10: (a) -\n{-\n[P(v x /a x )]} versus ln(v x /a x ) for N = 
350, <j) = 0.05 and rj = 0.9 (o), N — 500, = 0.05 and rj = 0.9 
(D), N = 350, = 0.05 and rj = 0.8 (o), N = 350, = 0.25 
and 7) = 0.9 (A). The solid lines correspond to the fit as made 
by Rouyer and Menon and has an exponent a — 1.52. The 
range of the solid lines corresponds to half the range used by 
Rouyer and Menon in their fit, but contains about 80% of 
their data points 
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observed range. When At is reduced a clear crossover 
develops. The behavior of the velocity distribution for 
velocities higher than the crossover velocity is more com- 
plicated than in boundary heating. There is also a sharp 
kink at the high-velocity end, that we have been unable 
to explain so far. 

This can explain the differences between velocity distri- 
butions for uniform and boundary heating. When heat- 
ing uniformly all the particles are heated every timestep 
At, while only a few collisions occur in that time. The 
system is dominated by heating and we see a velocity 
distribution with a constant exponent over the entire ob- 
served velocity range. If we are heating through a bound- 
ary, there are many collisions in the dense region in the 
center of the box, but only a few particles escape to the 
boundary to get heated. In this case, the system is dom- 
inated by dissipation and we see in general a velocity dis- 
tribution that strongly deviates from a Gaussian. This 
difference between uniform and boundary heating can be 
quantified by the ratio q = Nh/Nc, where Nh and Nq 
are the mean number of heatings and the number of col- 
lisions a particle has experienced, respectively. For a uni- 
formly heated gas of N = 350, 4> = 0-02 and ij = 0.4 we 
find that q = 120. For the same parameters but heat- 
ing through the boundary we find q = 0.08. For again 
the same parameters but using two-point heating with 
At = 1.0 we find q = 0.012. 

We tested the effect of changing q on velocity distribu- 
tions obtained with boundary heating. When increasing 
the number of particles TV or the area fraction <f>, the av- 
erage number of collisions increases. One can show in a 
mean field approximation that q ~ (TV</)) _1//2 . The aver- 
age distance a particle travels between collisions is given 
by Icoll ~ 1/0- The average distance a particles travels 
between heatings Iheat is given by the dimensions of the 
container. For a box of area A the average distance be- 
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FIG. 11: (a) P(v x /a x ) for different values of At. (b) 
— \n{—\n[P(v x /cr x )]} versus ln(v x /a x ). Data is taken for 
TV = 350, <t> = 0.02, r) = 0.4 and At = O.Ol(o), 0.03(D), 
0.05(o), 0.10(A), 0.30(<), 0.50(V) and 1.00(o). 



tween boundaries is given by Iheat ~ A 1 / 2 — (N/cj)) 1 / 2 . 
Finally, we know Njj/Nq ~ Icoll/lfieat- Our simulation 
obeys this approximation very well. In Fig. El we show 
velocity distributions for r\ = 0.8 and different combina- 
tions of TV and 4>. We measure the heating-dissipation 
ratio q in the simulation and show velocity distributions 
with the same q on top of each other. For q = 1.3 and 
q = 0.13 we find excellent collapse for different TV and 
<j>, even when we scale the system by a factor 8. For 
q = 0.013, where spatial correlations become very strong, 
we still find reasonable collapse. As we increase q we ob- 
serve the usual pattern, where a crossover appears in a 
distribution that was initially close to a Gaussian. This 
means that the velocity distributions are almost exclu- 
sively controlled by two parameters; the coefficient of 
restitution rj, which is a material parameters, and the 
heating-dissipation ratio q, which depends on experimen- 
tal conditions. 




ln(Vo x ) 

FIG. 12: Velocity distributions for different values of the 
heating-dissipation rate q. Distributions with the same q are 
shown on top of each other. (A) q — 1.3 and we show N = 100 
and <j> = 1 ■ 10" 3 (o), TV = 200 and <f> = 5 ■ 10 _4 (D), TV = 800 
and (f> = 1.25 ■ 10~ 4 (0). (B) q = 0.13 and we show TV = 100 
and 4> = 0.08(o), TV = 200 and <j> = 0.04(D), TV = 400 and 
<t> = 0.02(0). (C) q = 0.013 and we show TV = 100 and 
4> = 0.4(o), TV = 200 and (j> = 0.2(D), TV = 400 and <j> = 
0.1(<0). Inset: Heating-dissipation ratio q for TV = 800(o), 
TV = 400(D), TV = 200(<» and TV = 100(*) for several values 
of <f>. The line is a fit of the form (TV^) 1/2 . 



V. A SIMPLE MODEL WITHOUT SPATIAL 
DEGREES OF FREEDOM 

To study the dependence of the velocity distributions 
on the ratio q we constructed a model of the inelastic 
gas inspired by work done by Ulam 0] on elastic gases. 
Ulam rederived the Maxwell-Bolzmann distribution for 
the velocities of particles in a perfect gas using the fol- 
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lowing very simple model: he considered N particles with 
random initial velocities. At every time step a pair of 
particles was selected and the velocities of the particles 
were changed as if they had collided elastically. Ulam 
found that, regardless of the initial velocities, the system 
quickly evolved to a state where on average all energy was 
distributed equally over all particles and the deviation of 
the velocities of the particles around the average velocity 
of the system was given by the Maxwell-Bolzmann distri- 
bution. This approximation of random collisions is only 
justified when the gas is sufficiently dilute. The same 
approximation can also be made for the granular gas by 
replacing the elastic collisions between particles with in- 
elastic collisions. Since the gas cools down without any 
further energy input, heating has to be included as well. 
This allows us to study the dissipative behavior of the 
granular gas while neglecting all spatial correlations that 
develop in the more detailed simulations described be- 
fore. In particularly we can study the effect on the ve- 
locity distributions of changing the fraction q = Nh/Nq 
described in the previous section by varying the ratio of 
heatings and collisions. 

We adapted Ulam's procedure in order to model a two- 
dimensional granular gas consisting of N particles. Equa- 
tion Q] can be cast into the following form: 

, 1+7) ( cos 2 9 sin 9 cos 9 \ , . , 

V < = Vi - — { sinflcosfl sin 2 9 ) (v ' ~ ^ (4) 

where Vj and Vj are the velocities of particles i and j, rj 
is the coefficient of restitution and 9 is the angle between 
the separation vector and a reference axis. Collisions 
in our model occur by selecting at random particles i 
and j and an uniformly distributed impact parameter 
— 2R < b < 2R. We then use the above collision rule with 
9 — arcsin(6/2i?) + arccos(v- s/v), where v = (vj — Vj)/2 
is the velocity in the center-of-mass frame and s is a unit 
vector along the reference axis. We discard values of 9 
corresponding to (vj — v,) • r»j- < as these represent 
unphysical collisions. We heat the gas by selecting a 
random particle i and adding a random amount of ve- 
locity according to equation [21 To prevent the velocities 
from running away, we subtract the center-of-mass veloc- 
ity after heating. In a single time step, we have let 2C 
particles collide and we heat H particles. This gives us 
q = H/{2C). 

When we use inelastic instead of elastic collisions and 
drive the system, we find that the gas heats up until it 
reaches a steady state, where on average the energy lost 
in collisions is compensated by the heat inserted into the 
system through our driving mechanism. As H is cho- 
sen increasingly bigger in comparison to C, the granular 
temperature a of the gas increases accordingly. 

In figure IT31 we plotted the result for an inelastic gas 
with -q — 0.4. We varied the number of heatings and the 
number of collisions in a single time step from H = 100 
and C = 1 to H = 1 and C = 1000. As q is lowered, the 
velocity distributions develop a crossover and for g« 1 
the distributions are strongly non-Gaussian, similar to 



the velocity distributions obtained for At ^> 1 in two- 
point heating. In figure ITU we keep q — 0.025 fixed and 
vary 77. We see that the crossover point shifts down as 77 
is lowered and that the kink in the velocity distribution 
shifts down. 

We also compared velocity distributions found in the 
model with those acquired by simulation. To this end we 
measured q in simulations using different ways of heat- 
ing. For a uniformly heated gas of N = 350, <f> = 0.02 
and r\ = 0.4 we found q = 120. For the same parameters 
but heating through the boundary we found q = 0.08. 
For two-point heating with At — 1.0 we found q — 0.012. 
We ran the model with the corresponding set of parame- 
ters and measured velocity distributions. The results are 
shown in figure llSl 

The velocity distributions agree qualitatively, but espe- 
cially for the case of q=0.08 there is considerable quan- 
titative difference. Still, the model is able to generate 
non-Gaussian velocity distributions, while neglecting all 
spatial correlations. This suggests that the non-Gaussian 
behavior of the velocity distributions is not caused by 
spatial fluctuations and correlations in the gas. Instead 
it is the flow of energy through the system, mediated by 
the inelastic collisions, that determines the shape of the 
velocity distribution. This also means that there is no 
a priori difference between inhomogeneous and homoge- 
neous heating: the difference in the shape of the velocity 
distribution, particularly the occurrence of the crossover, 
is only a consequence of the fact that in inhomogeneous 
heating in general the number of collisions between par- 
ticles exceeds the number of particles being heated. 




FIG. 13: Velocity distribution for an inelastic gas of N = 500 
and 77 = 0.4. (a) P(v x /a x ). (b) —]n{—ln[P(v x /a x )]} versus 
]n(v x /crx). Data is for different values of q = 7^: q = 50 (o), 
9 = 5 (□), q = 1 (o), q = 0.5 (A), q = 0.05 (<), q = 5 • 10 -3 
(V) and q = 5 ■ 10~ 4 (>). 
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VI. RELATION TO EXPERIMENTS 

Due to the idealized nature of our system, it is not 
possible to do a direct comparison between our simula- 
tion and experiments. In the experiment of Rouyer and 
Menon the driving is through the boundary, but there are 
some significant differences between their heating mech- 
anism and the one we use in simulations with boundary 
heating. Due to gravity and the geometry of the setup 
the injection of energy in the experiment of Rouyer and 
Menon is mainly in the vertical direction. This energy 
is transferred into the horizontal direction by collisions 
between particles. Another difference is that in the ex- 
periment the frequency of driving is relatively low. Be- 
cause of this the dynamics of the gas close to the driving 
boundary is strongly dependend on the phase of the driv- 
ing cycle. In fact, it has been shown in simulation that 
for a system similar to the experiment by Rouyer and 
Menon, a Shockwave propagates up through the gas . 
At a certain distance from the boundary the time depen- 
dence has decayed away and the gas enters a steady state. 
It is in this steady state that the velocity distributions 
are measured. 

It is not yet established how this time dependence 
and the occurence of a shock wave influences the veloc- 
ity distributions in the steady state. A priori it is not 
clear if it is possible to compare velocity distributions in 
systems with a strong time-dependence, like the experi- 
ments, with those that have no time dependence, as is the 
case in our simulations. There are, however, reasons to 
assume this is possible. The velocity distributions in the 
experiment of Rouyer and Menon are measured only in 
the direction orthogonal to the driving direction. Simula- 
tions fijl show that the effect of the shock in the direction 




FIG. 14: Velocity distribution for an inelastic gas of N = 
500 and q = 0.025. (a) P(v x /a x ). (b) ~ln{-\n{P{v x /a x )]} 
versus \n(v x /a x ). Data is for r\ = 0.9 (o), 0.7 (□), 0.5 (o), 0.3 
(A) and 0.1 (<). 



orthogonal to the shock are usually relatively weak and 
decay rapidly in height. In the steady state, influence 
of the shock is absent in the orthogonal direction, even 
while it still may be apparent in the direction perpen- 
dicular to the driving direction. So, if we only look at 
velocity distributions in the orthogonal direction and in 
the steady state, a comparison between the experiment 
and our simulations is be valid. 

It is also not clear how in the experiments the dynam- 
ics of the gas shape the velocity distribution and whether 
it is controlled by the parameter q. We speculate that in 
the steady state the system behaves in fact like a one 
dimensional inelastic gas. Fast upward moving particles 
inject energy in the orthogonal direction when colliding 
with particles in the steady state, effectively function- 
ing as a heat source. In this picture the mean number 
of collisions between fast upwards moving particles and 
particles in the steady state would be Njj, the average 
number of heatings, and collisions between the particles 
in the steady state mutually would be Nc, the average 
number of collisions. One way of changing the shape of 
the velocity distribution would be changing the fraction 
of particles in the steady state. More particles in the 
steady state would lower Njj and increase Nc leading to 
more non-Gaussian velocity distributions. 

The above considerations not only apply to the exper- 
iment of Rouyer and Menon but to most of the other 
experiments as well. In the experiments of Blair and Ku- 
drolli and those of Olafsen and Urbach velocity distribu- 
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FIG. 15: Velocity distributions for an inelastic gas of N = 
350, 4> = 0.02 and r\ = 0.4. (a) P(v x /a x ) for q = 120 (solid), 
q = 0.08 (dashed) and q = 0.012 (dotted). Results on the 
bottom are obtained for simulation, results on the top for the 
model, (b) — ln{—\n[P(v x /a x )]} versus \n(v x /a x ). The sym- 
bols shown are velocity distributions acquired by simulation 
for q = 120 (o), 0.08 (□), 0.012 (o). The lines show the ve- 
locity distributions found in the model for the same values of 
q (solid, dotted, dashed). 
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tions are measured orthogonal to the driving direction. 
In both cases it is not so much the collisions with the 
bottom plate that drive the gas in the ortogonal direc- 
tions, but mainly collisions between fast upward moving 
particles with particles that have low velocities in the 
orthogonal directions. It is in these experiments rather 
than those of Rouyer and Menon that we find a similar 
dependence on 77 and q as we describe in this article. 

In the setup of Olafsen and Urbach [ll| velocity distri- 
butions go from non-Gaussian to Gaussian when a rough 
plate is used instead of a flat plate. On a flat plate, 
energy is only injected in the in-plane directions by off- 
angle collisions between neighboring particles. With a 
rough plate, energy is injected directly into the direc- 
tions parallel to the plate every time a particle collides 
with the plate, effectively increasing the number of heat- 
ings over collisions. Baxter and Olafsen |l2( observe the 
same behavior in a system where a layer of heavy parti- 
cles is inserted between the other layer of particles and a 
flat bottom plate. Particles from the upper layer have off- 
angle collisions with the layer of heavy particles, injecting 
energy in the in-plane directions every cycle. Particles 
in the upper layer show Gaussian velocity distributions, 
whereas particles in the lower layer have non-Gaussian 
velocity distributions. 

Most convincing is the experiment by Blair and Ku- 
drollli 8]. Here the number of collisions is increased by 
adding more particles. As a result, their velocity dis- 
tributions develop the same crossover we see both in our 
simulations and model. The reason why these transitions 
are not visible in the experiment by Rouyer and Menon 
as they increase number of particles, is that in the first 
case the effective coefficient of restitution is much lower, 
rj ~ 0.5, due to friction with the inclined plane. 

Again, one of the main problems considering the veloc- 
ity distributions in granular gases is that different setups 
and experiments usually find different velocity distribu- 
tions. As we have shown in this section, our finding of the 
controlling parameter q could ultimately explain these 
seemingly inconsistent results. 

VII. CONCLUSION 

We compared the velocity distributions of a granular 
gas that was driven by uniform heating and by heating 
through the boundary. Although it is usually assumed 
that uniform heating yields the same behavior as bound- 
ary heating when the boundary is far away and the spa- 
tial gradients are small, we find that there are clear qual- 
itative differences. When driven through the boundary, 
for instance, the gas can form coexisting cool liquid-like 
clusters surrounded by a hot gaseous state for certain 
values of <fi and rj. Such clusters do not occur in our 



simulations with uniform heating. 

The difference between uniform heating and boundary 
heating also extends to the velocity distributions. For 
uniform heating, we find velocity distributions that are 
close to Gaussian over the observed range of velocities 
and for a wide range of <f> and rj. When heating through 
the boundaries, on the other hand, we find that velocity 
distributions are stretched Gaussians that usually cross 
over from one exponent a to another for the high- velocity 
tail. For the high-velocity tail we find that the exponent 
a varies over a wide range. For strong dissipation this 
tail cannot be described by a single exponent. 

This means that there is no universal velocity dis- 
tribution with a = 1.5, as was proposed by Rouyer 
and Menon. Instead, the velocity distribution found by 
Rouyer and Menon may be just one out of the many 
distributions described here. The apparent universality 
may be due to the use of a narrow range of parameters 
(nearly elastic particles were used with rj = 0.93 and vary 
0.05 < <f> < 0.25). We find that varying cf> only has a mi- 
nor effect on the shape of the distribution when 77 is so 
close to the elastic limit. For low 77, we instead find ve- 
locity distributions that look much like the distributions 
found by Blair and Kudrolli for gases with 77 w 0.5. 

Furthermore, we demonstrate that the distribution of 
velocities for dissipative gases, while not universal in 
form, depends only on two parameters: the coefficient 
of restitution 77 (a material parameter) and q = Nh/Nq, 
the average ratio of heatings and collisions in the gas (a 
function of experimental conditions) . We find that veloc- 
ity distributions range from Gaussian for 5 > 1, where 
heating dominates dissipation, to strongly non-Gaussian 
for g « 1, where the dynamics of the gas is dominated 
by the dissipative collisions between particles. The dif- 
ferences in distributions we find between uniform and 
boundary heating can then be understood as different 
limits of q, for q 3> 1 and q < 1 respectively. We can 
control the parameter q more easily using two-point heat- 
ing and here we observe the transistion form Gaussian to 
strongly non-Gaussian behavior directly as we change q. 

Furthermore, a simple model of a driven, inelastic gas 
without spatial degrees of freedom reproduces the entire 
family of velocity distributions we find in simulation, as 
we vary 77 and q. This means that the velocity distri- 
butions are non-Gaussian not because of spatial corre- 
lations. Rather, it is the cascade of energy from a few 
high-energy particles to the slow-moving bulk of the gas 
that is the key determinant of the non-Gaussian veloc- 
ity distributions. These observations should aid in the 
construction of a kinetic theory of dissipative gases. 

We thank L. P. Kadanoff, N. Menon, A. Kudrolli, D. 
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